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Abstract 

Air showers recorded by the Haverah Park Array during the years 1974-1987 have been re¬ 
analysed. For the original estimate of the energy spectrum, a relationship between the ground 
parameter p(600) and the primary energy, as determined by Hillas in the 1970s, was used. Here 
we describe the energy spectrum obtained using the QGSJET98 interaction model in the CORSIKA 
Monte Carlo code, together with GEANT to simulate the detailed detector response to ground par¬ 
ticles. A new energy spectrum in the range 3 x 10^^ eV to 4 x 10^® eV is presented. 


1 Introduction 

After many years of research the origin of cosmic rays remains unclear. The composition of cosmic rays 
at high energies is not well known and even the energy spectrum is rather uncertain. Progress is limited 
primarily because cosmic rays above 10^® eV can only be observed, with reasonable statistics, through the 
extensive air showers of secondary particles that they produce in the atmosphere. The primary energy and 
mass of the initiating cosmic ray have to be deduced from the form and particle content of the showers. 
This, unfortunately, requires use of an air shower model, based on our knowledge of electromagnetic and 
hadronic particle interactions and particle decay and transport in the atmosphere. These models suffer 
from poor knowledge of the hadronic and nuclear interactions at high energies and so the results inferred 
from the measurements are, to some extent, model dependent. 

At Haverah Park, UK, a 12 km^ air shower array of water-Cherenkov detectors, was operational from 
1967-1987 to measure cosmic rays in the energy range 10^^ eV to 10^° eV. Here we present a re-analysis of 
data taken during 1974-1987. At that time the energy reconstruction of individual air shower events and 
the energy spectrum were obtained with a shower model of Hillas , that gave an empirical description 
of high-energy interactions. 

In recent years a variety of more sophisticated models for hadron production, based on Gribov-Regge 
theory, have become available which attempt to include quantum chromodynamics in a consistent way. 
The QJSJET98 model (Quark-Gluon-String model with jets) for high energy interactions, often used 
within the CORSIKA code developed for Monte Carlo calculations of air showers, has been rather 
successful in describing consistently a variety of experimental results over a wide energy range. 

Tools for detailed simulation of the response of the detectors to particles of different type, energy and 
impact angle have also become available, e.g. GEANT 0, that allow us to study details of the measuring 
process in a manner inconceivable 25 years ago. The aim of this paper is to present a re-analysis of 
Haverah Park data on the energy spectrum using the QGSJET98 model in the CORSIKA code and with 
GEANT simulations. In the following paper, in this volume, we present also new results on the cosmic 
ray mass composition 
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2 The Haverah Park Array 

The Haverah Park extensive air shower array was situated near Leeds, UK, at an altitude of 220 m above 
sea level, corresponding to a mean atmospheric depth of 1016 g cm“^, at 53° 58' N, 1°38' W. The particle 
detectors of this array were made up of water-Cherenkov modules of two types. The majority, which were 
used throughout the experiment, were galvanised iron tanks 2.29 m^ in area and filled to a depth of 1.2 
m with water brought from a nearby borehole in magnesium limestone. All tanks were lined with a white 
diffusing material (ICI Darvic). The water in each tank was viewed with an EMI photomultiplier (EMI 
9618YB) with a 5 inch diameter Sll Sb/CsO photocathode held so that it just dipped into the water. 
Large area detectors were achieved by grouping together a number of modules in huts with roofs having 
a thickness less than 4 g cm“^. The signals from the modules were added locally within each detector 
hut. 

The number of Cherenkov photons released in a water tank is closely proportional to the energy deposit 
of the shower particles in the water. Shower electrons and photons, with typical energies of 1-50 MeV, 
are effectively stopped, whilst muons (with a typical energy of 1-2 GeV) penetrate the tank and release 
a number of photons that is proportional to their track length in the water. As most of the energy of an 
air shower at ground is carried by the electromagnetic particles, the use of deep detectors is very effective 
for measuring the energy flow in the shower disc. The densities of Cherenkov photons per unit detector 
area (Cherenkov densities) were recorded in terms of the average signal from a vertical muon (1 vertical 
equivalent muon = 1 vem) per square metre. This signal corresponds to approximately 14 photoelectrons 
for Haverah Park tanks. Fig. shows the overall layout of the Haverah Park Array - a more detailed 
description of the array is given in . 

The signals from 15 of the 16 tanks in one A-site hut were summed to provide the signal used for triggering 
and for the density estimate. One tank in each hut was used to provide a low gain signal, used when the 
signal from the other 15 tanks is saturated. An overall trigger was formed if, in the central detector (Al) 
and in at least 2 out of the three remaining A-sites (A2-A4), a density of > 0.3 vem m“^ was recorded. 
The trigger rates were monitored daily over the life of the experiment. After correction for atmospheric 
pressure variations, the trigger rates were stable to better than 5%. 

In the energy range of interest here, the detectors located in a ring of 2 km from Al do not usually 
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Figure 1: Layout of the Haverah Park Array. A) The whole array. B) The orientations and relative heights of 
the detector huts A1-A4. C) The arrangement of water tanks within one of the four main A-site huts. 
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show signals above threshold. However they are important for the analysis because they constrain the 
air shower core position inside this ring. Relevant density information is also obtained from three 9 m^ 
detectors at 150 m from the center of the array. A relatively densely-packed sub array of 1 m^ detectors 
surrounding Al (not shown in Fig. |^) was not used in this re-analysis, as these detectors ran only for 4 
years. 

The experimental data for each event include the densities and arrival times of shower particles for 
each detector. The times were measured relative to the start of the pulse in Al, with an experimental 
measurement error of 20 ns. Typically in an event the 7 central detectors have signals recorded. For 
each event the arrival direction 0, 0, core position x,y, and the water-Cherenkov signal density at a core 
distance of 600 m, /9(600), were obtained. In total the database comprises 84455 events in the zenith 
angle range between 0° and 45° and with energies > 2 x 10^^ eV. For the present work these data were 
re-analysed using more detailed Monte Carlo calculations and improved reconstruction algorithms, as 
described below. 


3 Determination of Shower Direction, Core Position and Shower 
Size 

The Haverah Park data analysis procedure has three steps: firstly, the direction of the incoming cosmic 
ray is determined; secondly, the core location is estimated from the pattern of Cherenkov densities and 
/9(600) is determined; thirdly, the primary energy is estimated from the relationship between p(600) and 
the primary energy obtained by simulations. The zenith and azimuthal arrival directions are determined 
by fitting a plane to the particle arrival times recorded at the detectors Al to A4. The error in the zenith 
angle determination was estimated from the uncertainties introduced by the timing measurements. The 
absolute pointing accuracy was confirmed by comparisons with the directions of high energy muons 
observed in a magnetic spectrograph located at Al Q| and by comparison with estimates made with an 
air-Cherenkov array and with an array of scintillator detectors centered on Al Q. The zenith angle error 
is adequately represented by ag = 2.5° x sec0, for 0° < 0 < 45° and the azimuth angle is always smaller 
than that 1°). We use in this work the directions determined in the original Haverah Park analysis. 
Once the arrival direction has been found, the detector coordinates are projected into the shower plane, 
i.e. the plane normal to the shower axis passing through Al. Subsequently the core location (x,y) and 
/9(600) in the shower plane are obtained. For this purpose we have adopted a shape for the lateral 
distribution function that was directly measured 1^. The lateral distribution of the water-Cherenkov 
density p(r), in units of vem/m“^, in the core distance range r = 50-800 m is given by the function: 

p{r) = k r-(''+’'/4000 m) ^ ^ 

where r is in metres, fc is a normalization constant, and the slope parameter ry is given by 77 = 3.78 — 
1.44 sec 9. It has been shown that this dependence is in good agreement with simulation results. 
Due to the small number of densities usually available, it was not possible to fit values of ry reliably for 
each individual shower. However, the value of p(600) is only weakly dependent on ?y, which is one of 
the reasons for the choice of p(600) as the ground parameter from which we derive the primary energy. 
Specifically the mean value of p(600) increases (decreases) by only 15% if the reconstruction is made 
using a value of ry which is 0.3 larger (smaller) than the predicted value. Such a shift in the value of rj 
corresponds to ± 3 tr in terms of the intrinsic shower-to-shower fluctuations of ry ||^. We have ignored 
the dependence of ry on energy, which is measured to be 0.165 ± 0.022 per decade. 

The core position (x,y) and the shower size are fitted by comparing the measured densities with those 
predicted (from Eq. |^), i.e. by minimising: 

n 

= ^ — {k fin) - p{ri)f ( 2 ) 


3 


where pi are the measured densities at detector location i and Ui are uncertainties related to the density 
measurement and particle fluctuations. By differentiating Eq. the value of k which minimises is : 




( 3 ) 


The minimization technique is a special case of the maximum likelihood method of parameter esti¬ 
mation. If the probability distribution function for each observed density is a Gaussian with mean value 
Pi{r) and variance af then x^ is related to the likelihood L according to : 

X^ = —2 ln(L) -I- const. 

In this case, however, the weights depend on x, y, k and so Eq. ^ is no longer strictly valid. Another 
constraint on the formula is that it is applicable only if the measurement uncertainty is well described 
by a Gaussian distribution. When the number of expected particles is small it is more appropriate to 
use Poissonian statistics when calculating the likelihood function. Furthermore densities below trigger 
threshold or above saturation level cannot be used in the x^ method, but it is straightforward to incorpo¬ 
rate them in the likelihood method. As the measurement uncertainties are not known sufficiently well to 
justify a full maximum likelihood analysis we have adopted a compromise solution, as described below. 
For a given core position and zenith angle we assume a lateral distribution function (with the value of 77 
fixed to the mean value for that given zenith angle) and calculate the shower size constant k with Eq. 
We use only detectors with densities above threshold and below the saturation density. The variances ai 
used are obtained by adding in quadrature all the possible sources of error ||l^ : 
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where the first term is the absolute measurement uncertainty (Ej = (j{p)lp is the relative uncertainty), the 
second term is a constant uncertainty in p due to digitisation, and the third and fourth terms correspond 
to the Poisson uncertainty in the number of particles and Cherenkov photons, respectively. Ej and Xi are 
constants which depend on the type of the detector, Ai is the effective area of the detector, and pi is the 
“Poisson factor” - the effective number of real particles or photons that constitute the signal of a vertical 
equivalent muon. The Poisson factor can be > 1 because of the soft electromagnetic component and falls 
off with distance from the shower core because of the increasing number of muons relative to electrons 
and photons at greater distances. It is also dependent on zenith angle, as the muon to electromagnetic 
ratio varies with atmospheric depth. To investigate the uncertainty due to Poissonian fluctuations in 
the number of particles making up the water Gherenkov signal (as well as that due to fluctuations in 
the number of photoelectrons) signals from two adjacent detector modules were recorded for the same 
shower. The value of pi could be obtained with this investigation but limited statistics preclude an 
accurate determination of the dependence of pi on zenith angle and core distance. A general form for pi 
is given by: 

Pi = pIC + {I - p!C)N soft ( 5 ) 

where p/C is the fraction of the Cherenkov signal from muons, and iVgoft is the number of electromagnetic 
particles needed to make up the signal of a vertical muon. The value of pi used here was obtained from 
the parameterization of p/C as function of zenith angle and core distance in ||l^ with A^soft= 10 . 

Once we have obtained the value of fc, the predicted Cherenkov densities are calculated using pi = k /{r/) 
and converted to particle numbers by multiplying with pi. The likelihood function is calculated at 
this stage, including detectors above threshold and below saturation and using Poissonian or Gaussian 
statistics, depending on the particle number. A gradient search of the core position (x,y) is carried out to 
minimise the likelihood function. Once the optimum core position and k are found, the assumed lateral 
distribution function is modified by varying p around its expected value. For each value of 77 a gradient 
search is performed to obtain new values of k and core position. The whole procedure is repeated for 
varying starting values of the core position to avoid trapping in local minima. The minimum value of 
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Figure 2: Density map of a typical event recorded at the Haverah Park Array. The area of the circles indicates 
the size of the Cherenkov densities, which are shown also besides the detectors in vem m“^. Only the inner part 
of the array is displayed. The cross marks the derived core location. The p(600) of 0.55 vem m~^ at 19° of zenith 
corresponds to an energy of 0.3 EeV. 

the likelihood function for all values of rj, x, y, and k is assumed to be the true minimum. The value of 
p(600) is then obtained from the normalization constant k and the corresponding value of rj. 

Two main differences between this analysis and the original Haverah Park algorithm should be noted: pi 
varies with zenith angle and core distance and ry is allowed to vary in the fitting procedure, whereas both 
were considered constant originally. From the total event sample we re-analysed more than 8000 events 
with originally reconstructed core positions < 500 m away from Al. After the new core reconstruction 
a cut at 300 m was performed to ensure 100% trigger probability at energies above 3 x 10^^ eV for all 
particle types, allowing a simple and accurate calculation of the effective area for the reconstruction of 
the energy spectrum. 

Fig. H shows the density map of a typical event in the shower plane. Only the inner part of the array is 
displayed. 


4 Energy Calibration via Simulated Events 

4.1 Tank Response 

A detailed GEANT simulation of the propagation of electrons, gammas, and muons at different zenith 
angles through Haverah Park tanks has been performed. Cherenkov photons are ray-traced until they 
are absorbed or fall on the photocathode of the photomultiplier. The wavelength dependence of light 
absorption in the water and the reflectivity of the tank walls have been taken from fl^ but were normalised 
to the best estimates for the Haverah Park tanks. The peak values for Haverah Park are 83% for the 
wall reflectivity, 15 m absorption length in the water, and 22% photomultiplier quantum efficiency. The 
Thorn/EMI (9618YB PMT) wavelength dependent quantum efficiency was taken from the manufacturer’s 
specifications. The results of this simulation for particles that enter the tank vertically are summarised in 
Fig. The mean signal for a typical vertical muon « 1 GeV) was measured to be 14 photoelectrons 
[p^ , and this number is well reproduced by the simulations. 

By contrast, the mean energy of electrons and gamma rays in vertical showers are around 40 and 10 
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Figure 3: The mean signal produced in a Haverah Park tank by photons, electrons and mnons, that enter a tank 
vertically, as a fnnction of energy. Low energy electrons are absorbed in the lid of the tank. High energy muons 
can penetrate the tank and deposit « 240 MeV within it. 

MeV respectively at about 100-1000 m from the shower core. These calculated values are consistent 
with experimental evidence from Haverah Park work . Convolving the energy spectrum of electrons 
and gamma rays with the response given in Fig. yields mean signals of 2.6 and 0.9 photoelectrons for 
electrons and gamma-rays, respectively, at a typical core distance of 600 m. 

Specific simulations for different zenith angles (0°, 15°, 26°, 40° and 45°) and azimuth angles (0°, 45°, 90°) 
were performed to search for any change of the signal of electromagnetic particles due to border effects. 
For the mean photon signal, variations with angle were found to be less than 15%, and for electrons 
the variations are even smaller. The electromagnetic particles usually get completely absorbed in the 
tanks and the output signal is just proportional to the input particle energy. Thus, their contribution 
to the total signal at larger zenith angles is suppressed compared to muons because of the reduction of 
the projected area of the detectors. This effect can be included easily in the simulations if the mean 
signal from an electromagnetic particle at a given energy is independent of the arrival direction. The 
measured electromagnetic signal can also vary within a hut because of the varying azimuthal orientation 
of the tanks (see Fig. |^). For a given zenith angle the mean signal of the electromagnetic particles for 
the three azimuth angles was simulated. Using this mean value, the error in the electromagnetic signal 
when considering any azimuth is less than 7%. 

The signal due to muons, is proportional to the track length. For a given muon density the signal is 
also proportional to the tank area and as a result the mean signal is proportional to the tank volume 
and independent of the arrival direction of the shower relative to the tank. At large zenith angles the 
muon signal comprises a smaller number of muons (due to projected area) but with longer track length. 
Therefore Poisson fluctuations in the total number of muons going through a tank become more important 
at large zenith angles and this is accounted for in the simulations. 

To give an idea about the relative importance of the different contributions to the total water-Cherenkov 
signal, we have found from simulations that the ratio of the signal from muons to the total water- 
Cherenkov signal is 0.55 (0.42) at 600 m from the shower core for an iron (proton) shower arriving at a 
zenith angle of 26° with an energy of 0.4 EeV. 

In the following section we use the mean signals from muons, electrons and photons at different energies 
and angles to produce the water-Cherenkov lateral distribution function for individual events from Monte 
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Carlo simulations. The response of the tank to a particle at a given energy and arrival direction follows a 
distribution that is not always Gaussian. Although it would be more accurate to sample this distribution 
in the simulations to obtain the signal produced by a particle, we use the mean value of this distribution. 
We have investigated these differences in the density range of interest for this work: the systematic effect 
was found to be less than a 1%. 

4.2 Air Shower Simulations 

The Haverah Park data were originally analysed by comparison with shower simulations of Hillas and 
collaborators |^, using an ad-hoc hadronic interaction model. Subsequently accelerator measurements 
and the theoretical understanding of soft hadronic interactions have much improved, and the computing 
power available today allows much more detailed models to be used. At present the CORSIKA program 
is a standard tool for air shower simulations, and it is used by a variety of different experiments over 
the energy range from 1-10^^ GeV. CORSIKA uses EGS4 (Electron Gamma Shower code) for the 
simulation of the electromagnetic interactions, and features a detailed simulation of particle propagation 
and decay. Different models for the simulation of hadronic interactions are available within CORSIKA. 
For this work the QGSJET98 Q model was selected. QGSJET is based on Gribov-Regge theory of multi- 
Pomeron exchange to model high-energy soft hadronic interactions, at present the only viable theoretical 
approach to model cosmic ray interactions in the atmosphere. The QGSJET model also contains a 
treatment of hard collisions and the production of mini-jets, which become dominant at very high energy 
collisions, and allows a realistic simulation of nucleus-nucleus collisions. The interaction parameters are 
tuned to fit a wide range of accelerator results and to provide a secure extrapolation to higher energies. 
CORSIKA/QGSJET is fast and is able to describe a wide range of experimental cosmic ray results: from 
Cherenkov telescopes at 10^^ eV ||T^, over measurements of the hadronic, muonic and electromagnetic 
shower components in the knee region |l^, to lateral distributions and arrival times at « 5 x 10^^ eV Q, 
and air showers at 10^° eV |Q. 

Showers initiated by primary protons and iron nuclei were simulated at zenith angles of 0°, 15°, 26°, 40° 
and 45° with an energy of 4 x 10^^ eV. For a zenith angle of 26° showers with 0.2, 0.4, 0.8, 1.6, and 3.2 
EeV were also produced. Statistical thinning of shower particles was applied at the level of 10“®x Eq 
with a maximum particle weight limit of 10“^^ x Eo/eV. For each set 100 showers have been simulated, 
yielding a total of 3600 showers. 

The CORSIKA output was convolved with the tank response described in the previous section. We then 
fitted the resulting lateral distribution function for muons and electromagnetic particles to the function 
given in Eq. ^ for each individual shower. To increase the statistics further each shower was thrown 100 
times on to the array with random core positions ranging out to 400 m from the centre of the array. The 
zenith angles for these multiple Monte Carlo events is obtained by fluctuating the zenith angle of their 
parent event with the uncertainties described in the previous section. The muon and electromagnetic 
densities, including the Poissonian fluctuations, at the location of each of the detectors is obtained with 
the parametrizations described above. A check is made as to whether a simulated event fulfills the array 
trigger conditions and, if so, the densities are fluctuated according to measurement errors and recorded 
in the same format as real data. Each simulated event is then analysed with the algorithm described in 
the previous section to obtain the core position and p(600). Events with reconstructed cores at distances 
larger than 300 m from the central detector are rejected. It was found that the number of events in which 
the cores are falsely reconstructed to be outside 300 m is smaller than the number of events in which 
cores are falsely reconstructed to be inside 300 m. However, this is effect is small: at 0.4 EeV an increase 
of less than 1% in the number of events that pass the core distance cut is obtained, rising to 2% at 3.2 
EeV. 

4.2.1 Energy Calibration 

The correlation between /5(600) and the primary energy is used to determine the energy of an event. 
We have studied this correlation at a zenith angle of 26° since this is close to the median angle of the 


7 


I 

E 


o 

o 

to 




0.2 

0.18 

0.16 

0.14 

0.12 

0.1 

0.08 

0.06 

0.04 

0.02 


40 % 



Eo (EeV) 


Figure 4: Left panel: p(600) as a function of primary energy for showers intiated by proton and iron nuclei at 
6 = 26°, from CORSIKA/QGSJET simulations. The relation used for previous Haverah Park analyses is also 
plotted as a solid line. Right panel: the ratio between p(600) obtained from CORSIKA/QGSJET and by Hillas 
et al. Q 
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C 

1 

0.998±0.008 

0.616±0.006 

56 

1.011±0.004 

0.574±0.003 

Hillas 

1.018 

0.635 


Table 1: Parameters for relation between p(600) and primary energy for proton and Iron showers as calculated 
with CORSIKA/QGSJET. The parameters calculated by Hillas are also shown. Hillas used a reference zenith 
angle of 15°, so for comparison we have converted the parameters to a reference zenith angle of 26° using an 
attenuation length of 760 g cm“^. 

showers in the data. The value of p(600) for each energy was obtained by calculating the mean value 
of the reconstructed /9(600) for all simulated events generated for the corresponding set of CORSIKA 
showers. The results for proton and iron primaries are plotted in Fig. Q, and are compared with the 
result obtained by Hillas et al. [Q and used in previous Haverah Park analyses, e.g. |^. /9(600) grows 
nearly linearly with Ag, reconhrming that p(600) is a good energy indicator. The values of p(600) for 
iron-initiated showers are ^ 10 % higher than those for proton-initiated showers of the same energy, almost 
independent of energy. 

There are significant differences between the old and the new calibrations, which leads to « 30% lower 
reconstructed energies at about 10^® eV. The right panel of Fig. || shows the ratio between the old 
and two new calibrations, respectively. The slopes for the proton and iron calibrations are only slightly 
different, see Table ^ 

The relation between /9(600) and energy can be described by 

Fig (EeV) = Cp(600,26°)'>' ( 6 ) 

with C and 7 given in Table for proton and iron primaries, and for the Hillas model adopted in previous 
works. 
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Figure 5: Attenuation of p(600) with zenith from the constant intensity cut method, as obtained previously 
and in the present work. 

4.2.2 Attenuation Length 

Extensive air showers of different angles of incidence must be combined to derive an energy spectrum. 
Here, as the calibration has been computed for a zenith angle of 26°, the observed density, p{r,6), at 
zenith angle 9 and core distance r is corrected to that at 26° using the relation: 

p{r, 26°) = p{r, 9) exp 

where Aobs is the atmospheric depth at the observation site and A is the attenuation length with at¬ 
mospheric depth. A depends on the altitude of the array and on the type of detectors used. It can be 
determined from the experimental data and by simulation: A is usually obtained experimentally with the 
“constant intensity cut method”. This method is based on the assumption that the energy spectrum of 
cosmic rays is not dependent on the zenith angle of observation. p(600) spectra are obtained for a set 
of zenith angle bins. A horizontal slice of these spectra, at a fixed detection rate, then delineates the 
attenuation of p(600) with atmospheric depth. Fig. ^summarises the results from the constant intensity 
cut method on Haverah Park data as obtained by Edge et al. |^. Using a cut at an intensity of 10“^^ 
m“^ s“^ sr“^ (corresponding to an energy of about 10^® eV) in the zenith angle range 0° to 60° a value 
of A = 760 ± 40 g cm“^ was obtained. The slope of the line was calculated using a linear least square 
fit giving the points equal weight. We have repeated this analysis for the data accumulated between 
1974 and 1987, but restricted the analysis to zenith angles in the range 0°-45° and with the individual 
errors of the respective points taken into account in the fit. Our results on p(600) are compatible with 
the previous Haverah Park results at smaller zenith angles. Our estimate of the attenuation length is, 
however, 580 ± 50 g cm“^ and is shown in Fig. || as the dashed line. The disagreement comes mainly 
from the different zenith angle interval used in each calculation. The attenuation of p(600) with zenith 
angle is a result of the attenuation of the electromagnetic and muonic components; the muonic compo¬ 
nent attenuates more slowly, so when the contribution to /9(600) from the muonic component starts to 
dominate the attenuation length will increase. This becomes apparent at large zenith angles (see Fig. 
and therefore this analysis is restricted to zenith angles < 45°. 

The next step is to compare the attenuation length obtained from the data directly with the Monte Carlo 
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Figure 6: Attenuation of p(600) with zenith angle for different primaries and energies from Monte Carlo calcu¬ 
lations. The errors bars correspond to the spread of the value of p(600), to illustrate the shower fluctuations. 


Energy 

A 

A (g cm 

0.4 EeV 

1 

512±50 


56 

581±20 


Table 2: The attenuation length of p(600) determined for proton and iron primaries by simulation. 

simulations. For this purpose we have used the proton and iron-initiated showers at an energy of 0.4 EeV 
for different zenith angles. The results from simulating events to obtain mean values of the reconstructed 
p(600) are plotted in Fig. The attenuation lengths obtained from this figure are given in Table ||. 
They are consistent with the values from the re-analysis of the experimental data within the statistical 
uncertainties in the measurements and in the calculations. 

4.2.3 Energy Resolution 

The resolution of the energy reconstruction affects directly the energy spectrum obtained from data, 
shifting the spectrum upward if the energy resolution does not change with energy and shifting and 
changing the slope if it does Q. We have estimated this energy resolution from the simulated events for 
which we know the input core positions and energies exactly. 

As different zenith angles are combined to obtain the primary energy spectrum it is necessary to take into 
account changes of the energy resolution with zenith angle. Therefore, we have used simulated events 
for a fixed energy and primary mass at all zenith angles available and have weighted them according to 
the number of events in each zenith angle range. The reconstructed values of p(600) were converted to 
an equivalent p(600,26°) using the attenuation lengths given in Table Fig. ^ shows the results of this 
calculation. The spread in the values of p(600,26°) has different sources: (i) experimental reconstruction 
errors, (ii) intrinsic fluctuations in /9(600) due to shower development, and (hi) zenith angle errors which 
affect the calculation of p(600,26°) when using the attenuation length. 

The relative widths (t(p(600))/p( 600) of the distributions in Fig. ^ correspond directly to the energy 
resolution a{E)/E and are summarised in Table ^ The better resolution for iron primaries is attributed 
to the smaller intrinsic shower fluctuations. The effect of the energy resolution on the intensities reported 
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Figure 7: Resolution of p(600) integrated over all possible zenith angles. The attenuation lengths obtained from 
simultations corresponding to each primary mass have been used. 


Energy 

A 

cr(E)/E 

0.4 EeV 

1 

17% 


56 

12 % 

6.4 EeV 

1 

12 % 


56 

8 % 


Table 3: Energy resolution (t(E)/E for different primary energies and masses, 
for the energy spectrum is less than 2%, so no correction has been performed. 

4.2.4 Core Position 

Fig. 1^ shows the quality of the reconstruction of the core position for simulated proton events at two 
energies. The left panel shows the distribution of differences between the distances to detector A1 for the 
real and the reconstructed core positions. The right panel shows the distribution of distances between 
the real and reconstructed core position. The widths of both distributions are ~ 15 m. 

5 Energy Spectrum 

The showers recorded with the Haverah Park array, in the energy range of interest here, provide a 
sufficiently large database to allow for strict selection criteria without compromising statistical accuracy. 
All showers recorded between January 1974 and August 1987 were subject to the following conditions. 
Showers were required to have /o(600,26°) > 0.5 vem/m^. This cut removed all the low energy events 
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Figure 8: Right panel: Distribution of the differences between real (solid) and reconstructed (dashed) core 
position in m. Left panel: Distribution of the differences between the distances to the central triggering detector 
for real and reconstructed core position in m. 

where no reliable energy reconstruction was possible. It rejected about 90% of the event sample. Also 
shower cores had to land within 300 m of the array centre and zenith angles had to be less than 45°. 
These cuts reduced the sample by another 50%. The strict selection criteria are crucial, because they 
ensure 100% trigger efficiency over the full collection area of 300 m radius, independent of assumptions 
on the form of the lateral distribution function or primary particle type. The total on-time for these 
events was 3.657 x 10® seconds, about 11 years. The total aperture is 7.39 x 10^^ m^ s sr. 

A total of 3500 events survived all selection cuts and have been used to determine the energy spectrum. 
This number shows variations up to 2% for the different attenuation lengths obtained in simulations, and 
an increase of 30% if the attenuation length previously used for Haverah Park analysis (A = 760 ± 40 g 
cm“^) was adopted. 

The differential energy spectra obtained for different assumptions about the primary mass composition 
are listed in Table ^ and shown in Fig. In this graph we also show the energy spectrum obtained 
if the values of A obtained from the simulations are shifted by 40 g cm“^. This shift does not affect 
significantly (< 3%) the final energy spectrum, showing that the uncertainty of A has a negligible effect 
on our results. We use for the attenuation length the values of 512 g cm“^ for proton and 580 g cm“^ 
for iron to produce the final spectrum. 

Fig. |l^ shows the differential energy spectrum obtained for the assumption of pure proton and pure iron 
composition, respectively. The Haverah Park flux reported previously M is also shown. 

A least square fit to the experimental points in the energy range 3 x 10“ eV to 4 x 10^® eV is consistent 
with a power law of the form dN/dE = J (Eo/10^®elA)"’' with the index and normalization as given in 
Table |. 

It should also be noted that the difference in the if-p(600) relation between different primaries induces 
a difference of « 20% in the energy spectrum. In ref. ||^ it is shown that the lateral distribution data 
are in agreement with a bi-modal composition of 34% proton and 66% iron in the energy range 0.3-1.0 
EeV. Fig. 1^ shows the energy spectrum obtained if we adopt this composition for the cosmic ray beam. 
In this graph we also show a recent representation of the spectrum ||2^ derived from Akeno and Haverah 
Park data, and also the recent results obtained from the HIRes-MIA experiment . 


6 Conclusions 

The energy spectrum obtained in this analysis shows differences of up 30% at a given energy with a recent 
spectrum derived from Akeno and Haverah Park data. Our results are in good agreement with the 
recent results from the HIRes-MIA experiment ||^ in which a pure iron composition was assumed and 
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Figure 9: Energy spectra for different assumptions on primary mass and attenuation lengths as obtained from 
Haverah Park data. 



Figure 10: Energy spectrum obtained from Haverah Park for two assumptions of primary composition using 
CORSIKA/QGSJET compared with the spectra obtained previously. 
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Figure 11: Energy spectrum obtained from Haverah Park data assuming a bi-modal composition of 34% proton 
and 66% iron as obtained in |^. Other spectra are also shown for comparison. The continuous line shows the fit 
to the spectrum obtained in this work, see table ^ 

with results reported at the Hamburg conference using monocular HiRes dataJM. 

The spectral index after the knee at 3 x 10^® eV is estimated to be 3.1± 0.02 ||25|, so a further steepening 


(E) 

(EeV) 

JxlO^*^ 

(eV“^ m“^ sr“^ s“^) 


0.36 

78.5 ± 1.7 

2052 

0.53 

22.8 ± 0.8 

860 

0.76 

7.05 ± 0.36 

384 

1.10 

1.97 ± 0.16 

155 

1.54 

0.55 ± 0.07 

63 

2.29 

0.16 ± 0.03 

27 

3.49 

0.042 ± 0.01 

10 

4.52 

0.014 ± 0.006 

5 

7.07 

0.008 ± 0.004 

4 

8.60 

0.0028 ± 0.002 

2 


(E) 

(EeV) 

JxlO^^ 

(eV“^ m“^ sr“^ s“^) 

^showers 

0.34 

80.9 ± 1.8 

1986 

0.49 

22.9 ± 0.8 

815 

0.70 

7.31 ± 0.39 

378 

1.03 

2.00 ± 0.17 

150 

1.44 

0.52 ± 0.07 

57 

2.15 

0.177 ± 0.03 

28 

3.33 

0.039 ± 0.01 

9 

4.19 

0.015 ± 0.007 

5 

7.04 

0.0124 ± 0.005 

6 


Table 4: Differential fluxes from Haverah Park, assuming a pure proton (left) and a pure iron composition (right), 
respectively. 


A 

Index 

JxlO^^ at 10^« eV 
(eV“^ m“^ sr“^ s“^) 

1 

3.33±0.04 

2.66±0.09 

56 

3.34±0.04 

2.14±0.09 

Mixture 

3.33±0.04 

2.35±0.09 


Table 5: Parameterization of the energy spectra as power law for three assumptions of mass composition: pure 
proton, pure iron, and 34% proton 66% iron. See text for more explanations. 
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of the primary energy spectrum between the knee and 3 x 10^^ eV is needed to explain the spectral index 
derived in this work. It was claimed in |23 that a break in the spectrum occurs at an energy « 3 x 10^^ 
eV. Unfortunately this is the energy threshold of the spectrum described here. 

We have re-calculated the cosmic ray flux corresponding to the 4 highest energy events obtaining: 
JE^=8.27l4'o X 10^"^ eV^ sr“^ m“^ at an energy of 7.62 x 10^® eV. The average energy of these four 
events that were before above 10^° eV, is shifted a « 30% downwards to energies below 10^° eV. The 
most energetic event has an energy of 8.3 x 10^® eV. 

The Haverah Park array can be considered as an early prototype of the Auger Observatory which will 
employ water Cherenkov tanks of identical depth. The steps taken to derive an energy spectrum for 
the Auger Observatory are likely to be similar to those described here, and similar problems, like the 
attenuation length to be used and the ^-^(lOOO) relation, will need to be addressed. 
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